require(phyloseq)
require(tidyverse)
require(phyloseq)
require(reshape2)
require(dplyr)
require(ggplot2)
require(microbiome)
Loading required package: microbiome
microbiome R package (microbiome.github.com)
Copyright (C) 2011-2020 Leo Lahti,
Sudarshan Shetty et al. <microbiome.github.io>
Attaching package: ‘microbiome’
The following object is masked from ‘package:micrUBIfuns’:
prevalence
The following object is masked from ‘package:psych’:
alpha
The following object is masked from ‘package:vegan’:
diversity
The following object is masked from ‘package:ggplot2’:
alpha
The following object is masked from ‘package:base’:
transform
require(vegan)
Load data
ps_dmn <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/ITS/DMN_ests_ITS.Rdata")
sample_data(ps_dmn)$Herbicide <- factor(sample_data(ps_dmn)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_dmn)$herb_time<-paste(sample_data(ps_dmn)$Herbicide, sample_data(ps_dmn)$Time, sep = "_")
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_dmn)$Mode<-sample_data(ps_dmn)$Herbicide
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")
sample_data(ps_dmn)$Mode<- as.factor(values[match(sample_data(ps_dmn)$Mode, index)])
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
sample_data(ps_dmn)$Herbicide <- as.factor(values[match(sample_data(ps_dmn)$Herbicide, index)])
ps_rare <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/ITS/HerbPt1_rare_ITS.Rdata")
sample_data(ps_rare)$Herbicide <- factor(sample_data(ps_rare)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_rare)$herb_time<-paste(sample_data(ps_rare)$Herbicide, sample_data(ps_rare)$Time, sep = "_")
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")
sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
sample_data(ps_rare)$Herbicide <- as.factor(values[match(sample_data(ps_rare)$Herbicide, index)])
ps_trans <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/ITS/HerbPt1_hel_trans_ITS.Rdata")
sample_data(ps_trans)$Herbicide <- factor(sample_data(ps_trans)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_trans)$herb_time<-paste(sample_data(ps_trans)$Herbicide, sample_data(ps_trans)$Time, sep = "_")
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_trans)$Mode<-sample_data(ps_trans)$Herbicide
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")
sample_data(ps_trans)$Mode<- as.factor(values[match(sample_data(ps_trans)$Mode, index)])
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
sample_data(ps_trans)$Herbicide <- as.factor(values[match(sample_data(ps_trans)$Herbicide, index)])
Remove samples that are outlines or under sequenced.
ps_dmn <- subset_samples(ps_dmn, sample_names(ps_dmn) != "G009SG")
ps_dmn <- subset_samples(ps_dmn, sample_names(ps_dmn) != "G095SG")
ps_dmn <- subset_samples(ps_dmn, sample_names(ps_dmn) != "G123SG")
ps_dmn <- subset_samples(ps_dmn, sample_names(ps_dmn) != "G129SG")
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G009SG")
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G095SG")
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G123SG")
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G129SG")
ps_rare_sub<-prune_taxa(taxa_sums(ps_rare) > 2, ps_rare)
ordinations and adonis testing with three separate objects (i.e., dmn, rarefied, transformed). Rare taxa are removed from rarefied and transformed to successfully ordinate. At this point, the transformed data will not ordinate. This section is full dataset ordinations.
ord_dmn<-ordinate(physeq = ps_dmn, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit = 1000)
ord_rare<-ordinate(physeq = ps_rare_sub, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit = 1000)
ps_trans_sub<-prune_taxa(taxa_sums(ps_trans) > 0.01, ps_trans)
ord_transformed<-ordinate(physeq = ps_trans_sub, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit = 1000)
create alphadiversity tables
sample_sums(ps_rare)
G001SG G002SG G003SG G004SG G005SG G006SG G007SG G008SG G010SG G012SG G013SG G014SG G016SG G017SG G018SG G019SG G020SG
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
G021SG G022SG G024SG G025SG G026SG G027SG G028SG G029SG G030SG G031SG G032SG G033SG G034SG G035SG G036SG G037SG G038SG
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
G039SG G040SG G041SG G043SG G044SG G045SG G046SG G047SG G048SG G050SG G051SG G052SG G053SG G054SG G056SG G057SG G058SG
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
G059SG G060SG G062SG G063SG G064SG G065SG G066SG G067SG G068SG G069SG G070SG G071SG G072SG G073SG G074SG G075SG G076SG
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
G077SG G078SG G079SG G080SG G081SG G082SG G083SG G084SG G085SG G086SG G087SG G088SG G090SG G091SG G093SG G094SG G096SG
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
G097SG G098SG G100SG G101SG G102SG G103SG G104SG G105SG G106SG G107SG G108SG G109SG G111SG G112SG G113SG G114SG G115SG
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
G117SG G118SG G121SG G122SG G124SG G125SG G126SG G127SG G128SG G130SG G131SG G132SG G133SG G134SG G135SG G136SG G137SG
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
G139SG G140SG G141SG G142SG G143SG G144SG G145SG G146SG G147SG G148SG G149SG G150SG G151SG G153SG G154SG G155SG G156SG
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
G157SG G158SG G159SG G160SG G161SG G162SG G163SG G164SG G166SG G167SG G168SG G169SG G170SG G171SG G172SG G173SG G174SG
1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
G176SG G178SG G179SG G180SG
1000 1000 1000 1000
ps_rare
phyloseq-class experiment-level object
otu_table() OTU Table: [ 2033 taxa and 157 samples ]
sample_data() Sample Data: [ 157 samples by 53 sample variables ]
tax_table() Taxonomy Table: [ 2033 taxa by 7 taxonomic ranks ]
alpha_div <- estimate_richness(physeq = ps_rare, measures = c("Observed", "Shannon", "Chao1"))
#pull out metadata and concatonate with alpha diversity metrics
md<-data.frame(sample_data(ps_rare))
alpha_div_md <- rownames_to_column(alpha_div, "Barcode_ID_G") %>% full_join(md)
Joining, by = "Barcode_ID_G"
alpha_div_md$Herbicide <- factor(alpha_div_md$Herbicide, levels = c("Non-Treated", "Handweeded", "Atrazine-Mesotrione", "Dicamba", "Glyphosate"))
Shannon Div plots - no significant differences among herbicide treatments at any of the three time points
ggplot(data = alpha_div_md, aes(Herbicide, Shannon, color= Herbicide)) + facet_grid(. ~ Time) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1) )
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_Shannon.pdf")
Saving 7.29 x 4.51 in image
aov_t1<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T1",])
plot(aov_t1$residuals)
summary(aov_t1)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 0.389 0.09720 1.26 0.299
Residuals 48 3.704 0.07716
aov_t2<-aov(Shannon~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T2",])
plot(aov_t2$residuals)
summary(aov_t2)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 0.1234 0.03086 0.508 0.73
Residuals 46 2.7936 0.06073
aov_t3<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3",])
plot(aov_t3$residuals)
summary(aov_t3)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 0.5376 0.13440 2.551 0.051 .
Residuals 48 2.5287 0.05268
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(aov_t3, "Herbicide")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3", ])
$Herbicide
diff lwr upr p adj
Handweeded-Non-Treated 0.086844630 -0.19738244 0.3710717 0.9078666
Atrazine-Mesotrione-Non-Treated 0.129708505 -0.14882205 0.4082391 0.6803851
Dicamba-Non-Treated -0.002861275 -0.30174872 0.2960262 0.9999999
Glyphosate-Non-Treated 0.273363070 -0.01086400 0.5575901 0.0647237
Atrazine-Mesotrione-Handweeded 0.042863876 -0.22867317 0.3144009 0.9914461
Dicamba-Handweeded -0.089705904 -0.38208716 0.2026754 0.9066040
Glyphosate-Handweeded 0.186518440 -0.09085878 0.4638957 0.3282566
Dicamba-Atrazine-Mesotrione -0.132569780 -0.41941651 0.1542769 0.6864852
Glyphosate-Atrazine-Mesotrione 0.143654565 -0.12788248 0.4151916 0.5679081
Glyphosate-Dicamba 0.276224345 -0.01615691 0.5686056 0.0724032
Adonis testing of herbicide treatments by time point
ps_adonis<-function(physeq){
otu_tab<-data.frame(phyloseq::otu_table(physeq))
md_tab<-data.frame(phyloseq::sample_data(physeq))
if(taxa_are_rows(physeq)== T){
physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
print(anova(vegan::betadisper(physeq_dist, md_tab$Herbicide)))
vegan::adonis(physeq_dist ~ Herbicide * Time + Total_Weed_Veg , data = md_tab, permutations = 1000)
}
#ps_adonis(ps_rare_sub)
#remove one sample with no vegetation measurement.
ps_rare_sub_57<-subset_samples(ps_rare_sub, sample_names(ps_rare_sub) != "G065SG")
ps_adonis(ps_rare_sub_57)
#ps_adonis(ps_trans_sub)
ps_dmn_57<-subset_samples(ps_dmn, sample_names(ps_dmn) != "G065SG")
#ps_adonis(ps_dmn)
ps_adonis(ps_dmn_57)
Ordination plots DMN by time point
ord_t1_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Run 0 stress 0.1335559
Run 1 stress 0.1309245
... New best solution
... Procrustes: rmse 0.06976384 max resid 0.2060325
Run 2 stress 0.1315098
Run 3 stress 0.132755
Run 4 stress 0.134753
Run 5 stress 0.1308758
... New best solution
... Procrustes: rmse 0.0538814 max resid 0.2344595
Run 6 stress 0.1331698
Run 7 stress 0.1323904
Run 8 stress 0.1346675
Run 9 stress 0.1350022
Run 10 stress 0.1322358
Run 11 stress 0.1334733
Run 12 stress 0.1315248
Run 13 stress 0.1331407
Run 14 stress 0.1327605
Run 15 stress 0.1332143
Run 16 stress 0.1369199
Run 17 stress 0.1339829
Run 18 stress 0.1331858
Run 19 stress 0.1360529
Run 20 stress 0.1352334
Run 21 stress 0.1308759
... Procrustes: rmse 0.001713885 max resid 0.006848176
... Similar to previous best
*** Solution reached
T1_dmn<-ggordiplots::gg_ordiplot(ord = ord_t1_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_dmn_T1.pdf")
Saving 7.29 x 4.51 in image
ord_t2_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Run 0 stress 0.1325262
Run 1 stress 0.1325269
... Procrustes: rmse 0.0001074309 max resid 0.000408716
... Similar to previous best
Run 2 stress 0.1325262
... New best solution
... Procrustes: rmse 0.001015874 max resid 0.005488789
... Similar to previous best
Run 3 stress 0.1327962
... Procrustes: rmse 0.01038273 max resid 0.06195753
Run 4 stress 0.1338926
Run 5 stress 0.1392078
Run 6 stress 0.1338929
Run 7 stress 0.1325264
... Procrustes: rmse 0.001066291 max resid 0.005401677
... Similar to previous best
Run 8 stress 0.1325261
... New best solution
... Procrustes: rmse 0.0001868934 max resid 0.001066974
... Similar to previous best
Run 9 stress 0.1334238
Run 10 stress 0.1325259
... New best solution
... Procrustes: rmse 0.0007505591 max resid 0.004597761
... Similar to previous best
Run 11 stress 0.132925
... Procrustes: rmse 0.01264762 max resid 0.07728014
Run 12 stress 0.1438152
Run 13 stress 0.132526
... Procrustes: rmse 0.0003451055 max resid 0.001743723
... Similar to previous best
Run 14 stress 0.1325259
... Procrustes: rmse 0.0003370863 max resid 0.00167214
... Similar to previous best
Run 15 stress 0.1334866
Run 16 stress 0.1325266
... Procrustes: rmse 0.0009514394 max resid 0.005508805
... Similar to previous best
Run 17 stress 0.1325265
... Procrustes: rmse 0.0008955128 max resid 0.005161325
... Similar to previous best
Run 18 stress 0.1449845
Run 19 stress 0.1342762
Run 20 stress 0.1325263
... Procrustes: rmse 0.0003284059 max resid 0.001374251
... Similar to previous best
*** Solution reached
T2_dmn<-ggordiplots::gg_ordiplot(ord = ord_t2_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_dmn_T2.pdf")
Saving 7.29 x 4.51 in image
#this time point needs to be checked out. The ordination is not working properly.
ord_t3_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Run 0 stress 0.1451733
Run 1 stress 0.1459287
Run 2 stress 0.1510555
Run 3 stress 0.1501906
Run 4 stress 0.1503806
Run 5 stress 0.1459215
Run 6 stress 0.1528692
Run 7 stress 0.1451733
... Procrustes: rmse 0.00281782 max resid 0.01088482
Run 8 stress 0.1502075
Run 9 stress 0.1507344
Run 10 stress 0.1478567
Run 11 stress 0.1490906
Run 12 stress 0.1451724
... New best solution
... Procrustes: rmse 0.0002730257 max resid 0.001361369
... Similar to previous best
Run 13 stress 0.1507808
Run 14 stress 0.146477
Run 15 stress 0.1510089
Run 16 stress 0.1507665
Run 17 stress 0.1518463
Run 18 stress 0.1459425
Run 19 stress 0.1451714
... New best solution
... Procrustes: rmse 0.001928599 max resid 0.006972218
... Similar to previous best
Run 20 stress 0.1459215
*** Solution reached
T3_dmn<-ggordiplots::gg_ordiplot(ord = ord_t3_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_dmn_T3.pdf")
Saving 7.29 x 4.51 in image
Ordination plots on rarefied data by time point.
ord_t1_rare<-ordinate(physeq = subset_samples(ps_rare_sub, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.15407
Run 1 stress 0.1540701
... Procrustes: rmse 0.0004569647 max resid 0.001896947
... Similar to previous best
Run 2 stress 0.1560276
Run 3 stress 0.1540705
... Procrustes: rmse 0.00348766 max resid 0.02182018
Run 4 stress 0.1560277
Run 5 stress 0.158209
Run 6 stress 0.1586751
Run 7 stress 0.1545097
... Procrustes: rmse 0.01243307 max resid 0.04244993
Run 8 stress 0.1540702
... Procrustes: rmse 0.003479511 max resid 0.02173529
Run 9 stress 0.1543042
... Procrustes: rmse 0.01763231 max resid 0.08581444
Run 10 stress 0.1541252
... Procrustes: rmse 0.007485728 max resid 0.03683984
Run 11 stress 0.1540706
... Procrustes: rmse 0.003606795 max resid 0.02250138
Run 12 stress 0.1548484
Run 13 stress 0.1541275
... Procrustes: rmse 0.00762305 max resid 0.03699455
Run 14 stress 0.1543157
... Procrustes: rmse 0.01699931 max resid 0.07997088
Run 15 stress 0.1560286
Run 16 stress 0.1560286
Run 17 stress 0.1581542
Run 18 stress 0.1540702
... Procrustes: rmse 0.003594385 max resid 0.02243185
Run 19 stress 0.1550319
Run 20 stress 0.1560267
*** Best solution repeated 1 times
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare$plot + theme_classic()
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare_plot<-T1_rare$plot + theme_classic() + xlim(-0.5, 0.5) + ylim(-0.5, 0.5) + guides(color=guide_legend("Treatment")) + xlab("NMDS 1") + ylab("NMDS 2")
T1_rare_plot
Warning: Removed 3 rows containing missing values (`geom_point()`).
library(cowplot)
my_legend <- get_legend(T1_rare_plot)
Warning: Removed 3 rows containing missing values (`geom_point()`).
library(ggpubr)
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordinationlegend.pdf")
Saving 7.29 x 4.51 in image
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordinationlegend.eps")
Saving 7.29 x 4.51 in image
T1_rare_plot<-T1_rare_plot + theme(legend.position = "none")
T1_rare_plot
Warning: Removed 3 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T1.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 3 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T1.eps")
Saving 7.29 x 4.51 in image
Warning: Removed 3 rows containing missing values (`geom_point()`).
ord_t2_rare<-ordinate(physeq = subset_samples(ps_rare_sub, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1716481
Run 1 stress 0.1720658
... Procrustes: rmse 0.04180674 max resid 0.1226574
Run 2 stress 0.1720724
... Procrustes: rmse 0.0407305 max resid 0.1223211
Run 3 stress 0.1701444
... New best solution
... Procrustes: rmse 0.08690941 max resid 0.2978527
Run 4 stress 0.1709893
Run 5 stress 0.1711196
Run 6 stress 0.1728237
Run 7 stress 0.1708441
Run 8 stress 0.1706909
Run 9 stress 0.1706258
... Procrustes: rmse 0.06561966 max resid 0.2792094
Run 10 stress 0.171183
Run 11 stress 0.1746478
Run 12 stress 0.1746525
Run 13 stress 0.1709914
Run 14 stress 0.1734335
Run 15 stress 0.1710644
Run 16 stress 0.1710837
Run 17 stress 0.1702885
... Procrustes: rmse 0.05689984 max resid 0.2828329
Run 18 stress 0.1734279
Run 19 stress 0.1715313
Run 20 stress 0.1713459
Run 21 stress 0.1729146
Run 22 stress 0.1736516
Run 23 stress 0.1703011
... Procrustes: rmse 0.05756458 max resid 0.2826533
Run 24 stress 0.170112
... New best solution
... Procrustes: rmse 0.005398775 max resid 0.02579526
Run 25 stress 0.1702884
... Procrustes: rmse 0.05862291 max resid 0.2821017
Run 26 stress 0.1706282
Run 27 stress 0.172058
Run 28 stress 0.174658
Run 29 stress 0.1706231
Run 30 stress 0.1702879
... Procrustes: rmse 0.05856164 max resid 0.2821153
Run 31 stress 0.1701561
... Procrustes: rmse 0.005990811 max resid 0.02698447
Run 32 stress 0.1715783
Run 33 stress 0.1757041
Run 34 stress 0.1709007
Run 35 stress 0.1715863
Run 36 stress 0.1742448
Run 37 stress 0.1702882
... Procrustes: rmse 0.05868747 max resid 0.2820668
Run 38 stress 0.1702889
... Procrustes: rmse 0.05873721 max resid 0.2820447
Run 39 stress 0.1715866
Run 40 stress 0.170299
... Procrustes: rmse 0.05905885 max resid 0.2819334
Run 41 stress 0.1709006
Run 42 stress 0.1702394
... Procrustes: rmse 0.01169203 max resid 0.03765
Run 43 stress 0.1767481
Run 44 stress 0.1720567
Run 45 stress 0.1711196
Run 46 stress 0.1703014
... Procrustes: rmse 0.05919363 max resid 0.2818825
Run 47 stress 0.1709005
Run 48 stress 0.1703017
... Procrustes: rmse 0.05921868 max resid 0.2818621
Run 49 stress 0.1709909
Run 50 stress 0.1745794
Run 51 stress 0.1709004
Run 52 stress 0.1710641
Run 53 stress 0.1731
Run 54 stress 0.1728168
Run 55 stress 0.1706225
Run 56 stress 0.1706709
Run 57 stress 0.1732789
Run 58 stress 0.1711185
Run 59 stress 0.1702091
... Procrustes: rmse 0.009167803 max resid 0.02910266
Run 60 stress 0.1703042
... Procrustes: rmse 0.05940232 max resid 0.2817527
Run 61 stress 0.1701124
... Procrustes: rmse 0.0009526391 max resid 0.002988642
... Similar to previous best
*** Best solution repeated 1 times
T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_rare_plot<-T2_rare$plot + theme_classic() + xlim(-0.5, 0.5) + ylim(-0.5, 0.5) + theme(legend.position = "none") + xlab("NMDS 1") + ylab("NMDS 2")
T2_rare_plot
Warning: Removed 4 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T2.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 4 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T2.eps")
Saving 7.29 x 4.51 in image
Warning: Removed 4 rows containing missing values (`geom_point()`).
ord_t3_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1685562
Run 1 stress 0.173731
Run 2 stress 0.1770624
Run 3 stress 0.1745126
Run 4 stress 0.1722637
Run 5 stress 0.1686696
... Procrustes: rmse 0.01027886 max resid 0.05714972
Run 6 stress 0.1722627
Run 7 stress 0.1684508
... New best solution
... Procrustes: rmse 0.01175199 max resid 0.06578898
Run 8 stress 0.1786699
Run 9 stress 0.1686779
... Procrustes: rmse 0.02156585 max resid 0.130108
Run 10 stress 0.1685888
... Procrustes: rmse 0.01036983 max resid 0.05925135
Run 11 stress 0.1769021
Run 12 stress 0.1698586
Run 13 stress 0.1798158
Run 14 stress 0.1774968
Run 15 stress 0.1723331
Run 16 stress 0.1734184
Run 17 stress 0.1704177
Run 18 stress 0.1686749
... Procrustes: rmse 0.02108582 max resid 0.1271145
Run 19 stress 0.1686812
... Procrustes: rmse 0.02201686 max resid 0.1329718
Run 20 stress 0.1684506
... New best solution
... Procrustes: rmse 0.0001109389 max resid 0.0005559503
... Similar to previous best
*** Best solution repeated 1 times
T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_rare_plot<-T3_rare$plot + theme_classic() + xlim(-0.5, 0.5) + ylim(-0.5, 0.5) + theme(legend.position = "none") + xlab("NMDS 1") + ylab("NMDS 2")
T3_rare_plot
Warning: Removed 4 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T3.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 4 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T3.eps")
Saving 7.29 x 4.51 in image
Warning: Removed 4 rows containing missing values (`geom_point()`).
library(ggpubr)
ggarrange(T1_rare_plot, T2_rare_plot, T3_rare_plot, ncol = 1)
Warning: Removed 3 rows containing missing values (`geom_point()`).
Warning: Removed 4 rows containing missing values (`geom_point()`).
Warning: Removed 4 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_ordination.pdf", width = 5, height = 10)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_ordination.eps", width = 5, height = 10)
CAP ordination plots rarefied
t1_dist <- distance(subset_samples(ps_rare, Time=="T1"), method="bray") #get wUnifrac and save
t1_table<-as.matrix(dist(t1_dist)) #transform wUnifrac index
ord_t1_rare_cap <- capscale(t1_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T1"))))
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T1_cap.pdf")
t2_dist <- distance(subset_samples(ps_rare, Time=="T2"), method="bray") #get wUnifrac and save
t2_table<-as.matrix(dist(t2_dist)) #transform wUnifrac index
ord_t2_rare_cap <- capscale(t2_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T2"))))
T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T2_cap.pdf")
#G166SG identified as outlier based on plots with it included. Removed to create plot.
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G166SG")
t3_dist <- distance(subset_samples(ps_rare, Time=="T3"), method="bray") #get wUnifrac and save
t3_table<-as.matrix(dist(t3_dist)) #transform wUnifrac index
ord_t3_rare_cap <- capscale(t3_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T3"))))
T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T3_cap.pdf")
Pairwise adonis testing no needed becasue of insignificant gloabal test.
ps_pairwiseadonis<-function(physeq){
otu_tab<-data.frame(phyloseq::otu_table(physeq))
md_tab<-data.frame(phyloseq::sample_data(physeq))
if(taxa_are_rows(physeq)== T){
physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
pairwiseAdonis::pairwise.adonis(x = physeq_dist, factors = md_tab$Herbicide, p.adjust.m = "none", perm = 1000)
}
ps_t1<-subset_samples(ps_rare_sub, Time == "T1")
ps_t1<-prune_taxa(taxa_sums(ps_t1) > 2, ps_t1)
ps_t2<-subset_samples(ps_rare_sub, Time == "T2")
ps_t2<-prune_taxa(taxa_sums(ps_t2) > 2, ps_t2)
ps_t3<-subset_samples(ps_rare_sub, Time == "T3")
ps_t3<-prune_taxa(taxa_sums(ps_t3) > 2, ps_t3)
ps_pairwiseadonis(ps_t1)
ps_pairwiseadonis(ps_t2)
ps_pairwiseadonis(ps_t3)
Pairwise betadispr by treatment, time and mode
ps_betadispr<-function(physeq, groupingvar = "Groupingvar"){
otu_tab<-data.frame(phyloseq::otu_table(physeq))
md_tab<-data.frame(phyloseq::sample_data(physeq))
if(taxa_are_rows(physeq)== T){
physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
mod<-vegan::betadisper(physeq_dist, md_tab[,groupingvar])
## Perform test
print(anova(mod))
## Permutation test for F
pmod <- vegan::permutest(mod, permutations = 1000, pairwise = TRUE)
print(pmod)
print(boxplot(mod))
}
permute test of dispersion
ps_betadispr(subset_samples(ps_rare_sub, Time == "T3"), groupingvar = "Herbicide")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 4 0.00635 0.0015867 0.1932 0.9408
Residuals 48 0.39424 0.0082133
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 4 0.00635 0.0015867 0.1932 1000 0.9451
Residuals 48 0.39424 0.0082133
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Atrazine-Mesotrione Dicamba Glyphosate Handweeded Non-Treated
Atrazine-Mesotrione 0.93007 0.53946 0.86913 0.7502
Dicamba 0.92524 0.42557 0.73726 0.5854
Glyphosate 0.54158 0.42422 0.62537 0.7832
Handweeded 0.84582 0.73548 0.61029 0.8442
Non-Treated 0.71761 0.59620 0.77349 0.82944
$stats
[,1] [,2] [,3] [,4] [,5]
[1,] 0.2201596 0.2491324 0.1949060 0.2643697 0.2208840
[2,] 0.2635563 0.2905369 0.2531917 0.2745103 0.2571909
[3,] 0.2896573 0.3043022 0.2848346 0.3109094 0.2923035
[4,] 0.3593706 0.3524222 0.3591617 0.3254448 0.3874820
[5,] 0.4966550 0.3757590 0.4464053 0.3722896 0.4658206
$n
[1] 12 9 11 11 10
$conf
[,1] [,2] [,3] [,4] [,5]
[1,] 0.2459558 0.2717092 0.2343518 0.2866449 0.2272048
[2,] 0.3333588 0.3368951 0.3353175 0.3351740 0.3574021
$out
G141SG G146SG G144SG
0.6092207 0.5282037 0.5427749
$group
[1] 1 2 4
$names
[1] "Atrazine-Mesotrione" "Dicamba" "Glyphosate" "Handweeded"
[5] "Non-Treated"
beta_boxplot<-function (physeq, method = "bray", group)
{
require("phyloseq")
require("ggplot2")
group2samp <- list()
group_list <- get_variable(sample_data(physeq), group)
for (groups in levels(group_list)) {
target_group <- which(group_list == groups)
group2samp[[groups]] <- sample_names(physeq)[target_group]
}
beta_div_dist <- phyloseq::distance(physeq = physeq, method = method)
beta_div_dist <- as(beta_div_dist, "matrix")
dist_df <- data.frame()
counter <- 1
for (groups in names(group2samp)) {
sub_dist <- beta_div_dist[group2samp[[groups]], group2samp[[groups]]]
no_samp_col <- ncol(sub_dist)
no_samp_row <- nrow(sub_dist)
for (cols in seq(no_samp_col)) {
if (cols > 1) {
for (rows in seq((cols - 1))) {
dist_df[counter, "sample_pair"] <- paste0(colnames(sub_dist)[cols],
"-", rownames(sub_dist)[rows])
dist_df[counter, "group"] <- groups
dist_df[counter, "beta_div_method"] <- method
dist_df[counter, "beta_div_value"] <- sub_dist[rows,
cols]
counter = counter + 1
}
}
}
}
plot_boxplot <- ggplot(data = dist_df, aes(x = group, y = beta_div_value,
color = group)) + geom_boxplot(outlier.shape = NA) +
geom_jitter() + theme_bw() + xlab(group) + ylab(method) +
theme(axis.text.x = element_text(angle = 45, vjust = 1,
hjust = 1))
list_Out <- list(data = dist_df, plot = plot_boxplot)
return(list_Out)
}
Box and whisker plots of distance within group distances
#remotes::install_github("antonioggsousa/micrUBIfuns")
#library(micrUBIfuns)
T1_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Herbicide")
T1_beta_plot <- T1_beta$plot
T1_beta_plot <- T1_beta_plot + theme_classic()+ guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75)
T1_beta_plot
Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).
my_legend <- get_legend(T1_beta_plot)
Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_beta_legend.pdf")
Saving 7.29 x 4.51 in image
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_beta_legend.eps")
Saving 7.29 x 4.51 in image
T1_beta_plot<-T1_beta_plot+ theme(legend.position = "none")
T1_beta_plot
Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T1_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T1_rare_withingroup_beta.eps")
Saving 7.29 x 4.51 in image
Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).
T1_beta_df<- T1_beta$data
T1_betamod<-aov(formula = beta_div_value ~ group ,data = T1_beta_df)
summary(T1_betamod)
Df Sum Sq Mean Sq F value Pr(>F)
group 4 0.2967 0.07419 8.83 1.1e-06 ***
Residuals 250 2.1003 0.00840
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(x = T1_betamod, which = "group")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ group, data = T1_beta_df)
$group
diff lwr upr p adj
Dicamba-Atrazine-Mesotrione -0.027993939 -0.07861966 0.022631781 0.5508308
Glyphosate-Atrazine-Mesotrione -0.036412121 -0.08703784 0.014213599 0.2806138
Handweeded-Atrazine-Mesotrione -0.058666667 -0.11176337 -0.005569963 0.0220819
Non-Treated-Atrazine-Mesotrione 0.040006061 -0.01061966 0.090631781 0.1940292
Glyphosate-Dicamba -0.008418182 -0.05644596 0.039609594 0.9889807
Handweeded-Dicamba -0.030672727 -0.08129845 0.019952993 0.4576855
Non-Treated-Dicamba 0.068000000 0.01997222 0.116027775 0.0012035
Handweeded-Glyphosate -0.022254545 -0.07288027 0.028371175 0.7468226
Non-Treated-Glyphosate 0.076418182 0.02839041 0.124445957 0.0001747
Non-Treated-Handweeded 0.098672727 0.04804701 0.149298448 0.0000019
T2_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Herbicide")
T2_beta_plot <- T2_beta$plot
T2_beta_plot <- T2_beta_plot+ theme_classic() + theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75)
T2_beta_plot
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T2_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T2_rare_withingroup_beta.eps")
Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).
T2_beta_df<- T2_beta$data
T2_betamod<-aov(formula = beta_div_value ~ group ,data = T2_beta_df)
summary(T2_betamod)
Df Sum Sq Mean Sq F value Pr(>F)
group 4 0.080 0.019995 4.367 0.00201 **
Residuals 231 1.058 0.004578
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(x = T2_betamod, which = "group")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ group, data = T2_beta_df)
$group
diff lwr upr p adj
Dicamba-Atrazine-Mesotrione -0.029188889 -0.070787897 0.0124101197 0.3048204
Glyphosate-Atrazine-Mesotrione -0.020375758 -0.057770485 0.0170189703 0.5647456
Handweeded-Atrazine-Mesotrione 0.022206061 -0.015188667 0.0596007884 0.4780643
Non-Treated-Atrazine-Mesotrione -0.015555556 -0.054775477 0.0236643659 0.8113359
Glyphosate-Dicamba 0.008813131 -0.031069709 0.0486959716 0.9738007
Handweeded-Dicamba 0.051394949 0.011512109 0.0912777898 0.0043255
Non-Treated-Dicamba 0.013633333 -0.027965675 0.0552323419 0.8962631
Handweeded-Glyphosate 0.042581818 0.007106064 0.0780575719 0.0097763
Non-Treated-Glyphosate 0.004820202 -0.032574526 0.0422149299 0.9966025
Non-Treated-Handweeded -0.037761616 -0.075156344 -0.0003668883 0.0465006
T3_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Herbicide")
T3_beta$plot #+ scale_color_manual(values = c("#F8766D", "#A3A500", "#00BF7D", "#00B0F6", "#E76BF3")) +
T3_beta_plot <- T3_beta$plot
T3_beta_plot <- T3_beta_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("")
T3_beta_plot <-T3_beta_plot + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T3_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T3_rare_withingroup_beta.eps")
Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
T3_beta_df<- T3_beta$data
T3_betamod<-aov(formula = beta_div_value ~ group ,data = T3_beta_df)
summary(T3_betamod)
Df Sum Sq Mean Sq F value Pr(>F)
group 4 0.0672 0.016792 1.689 0.153
Residuals 252 2.5054 0.009942
TukeyHSD(x = T3_betamod, which = "group")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ group, data = T3_beta_df)
$group
diff lwr upr p adj
Dicamba-Atrazine-Mesotrione 0.009452020 -0.04731207 0.06621612 0.9909320
Glyphosate-Atrazine-Mesotrione -0.037518182 -0.08753732 0.01250096 0.2404830
Handweeded-Atrazine-Mesotrione -0.014736364 -0.06475550 0.03528278 0.9275528
Non-Treated-Atrazine-Mesotrione -0.023409091 -0.07637301 0.02955483 0.7430333
Glyphosate-Dicamba -0.046970202 -0.10570358 0.01176317 0.1840580
Handweeded-Dicamba -0.024188384 -0.08292176 0.03454499 0.7896781
Non-Treated-Dicamba -0.032861111 -0.09412180 0.02839957 0.5804264
Handweeded-Glyphosate 0.022781818 -0.02946147 0.07502511 0.7524585
Non-Treated-Glyphosate 0.014109091 -0.04096017 0.06917835 0.9554583
Non-Treated-Handweeded -0.008672727 -0.06374199 0.04639653 0.9926689
library(ggpubr)
ggarrange(T1_beta_plot, T2_beta_plot, T3_beta_plot, ncol = 1)
Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_within_group_beta.pdf", width = 7, height = 10)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_within_group_beta.eps", width = 7, height = 10)
Examination of dissimilarity across time points by treatment and then again by all chemical treatments combined.
T1_beta_df$Time<-"T1"
T2_beta_df$Time<-"T2"
T3_beta_df$Time<-"T3"
beta_div_T1_T2_T3 <- rbind(T1_beta_df, T2_beta_df, T3_beta_df)
beta_anova<-function(data, Herbicide = "Herbicide"){
df_sub<- data %>% filter(group == Herbicide)
mod<-aov(beta_div_value ~ Time, data = df_sub)
print(summary(mod))
print(TukeyHSD(mod, "Time"))
boxplot(df_sub$beta_div_value ~ df_sub$Time)
}
#beta_anova(beta_div_T1_T2_T3, Herbicide = "Non-Treated")
#beta_anova(beta_div_T1_T2_T3, Herbicide = "Handweeded")
#beta_anova(beta_div_T1_T2_T3, Herbicide = "Dicamba")
#beta_anova(beta_div_T1_T2_T3, Herbicide = "Atrazine-Mesotrione")
#beta_anova(beta_div_T1_T2_T3, Herbicide = "Glyphosate")
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide
index <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Handweeded", "Non-Treated")
sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])
#+ scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T1_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Mode")
T1_beta_chemical_combined_plot <- T1_beta_chemical_combined$plot
T1_beta_chemical_combined_plot<- T1_beta_chemical_combined_plot + theme_classic() + guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T1_beta_chemical_combined_plot
Warning: Removed 9 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 9 rows containing missing values (`geom_point()`).
my_legend <- get_legend(T1_beta_chemical_combined_plot)
Warning: Removed 9 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 10 rows containing missing values (`geom_point()`).
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_beta_combined_legend.pdf")
Saving 7.29 x 4.51 in image
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_beta_combined_legend.eps")
Saving 7.29 x 4.51 in image
T1_beta_chemical_combined_plot<-T1_beta_chemical_combined_plot+ theme(legend.position = "none")
T1_beta_chemical_combined_plot
Warning: Removed 9 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 10 rows containing missing values (`geom_point()`).
T2_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Mode")
T2_beta_chemical_combined_plot <- T2_beta_chemical_combined$plot
T2_beta_chemical_combined_plot<- T2_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T2_beta_chemical_combined_plot
Warning: Removed 3 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 3 rows containing missing values (`geom_point()`).
T3_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Mode")
T3_beta_chemical_combined_plot <- T3_beta_chemical_combined$plot
T3_beta_chemical_combined_plot<- T3_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T3_beta_chemical_combined_plot
Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).
ggarrange(T1_beta_chemical_combined_plot, T2_beta_chemical_combined_plot, T3_beta_chemical_combined_plot, ncol = 1)
Warning: Removed 9 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 9 rows containing missing values (`geom_point()`).
Warning: Removed 3 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 4 rows containing missing values (`geom_point()`).
Warning: Removed 5 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 5 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_within_group_beta_chemical_combined.pdf", width = 5, height = 10)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_within_group_beta_chemical_combined.eps", width = 5, height = 10)
T1_beta_df_chemical_combined <- T1_beta_chemical_combined$data
T2_beta_df_chemical_combined<- T2_beta_chemical_combined$data
T3_beta_df_chemical_combined<- T3_beta_chemical_combined$data
T1_beta_df_chemical_combined$Time<-"T1"
T2_beta_df_chemical_combined$Time<-"T2"
T3_beta_df_chemical_combined$Time<-"T3"
m1<-aov(beta_div_value ~ group, data = T1_beta_df_chemical_combined)
summary(m1)
TukeyHSD(m1, "group")
boxplot(beta_div_value ~ group, data = T1_beta_df_chemical_combined)
m2<-aov(beta_div_value ~ group, data = T2_beta_df_chemical_combined)
summary(m2)
TukeyHSD(m2, "group")
boxplot(beta_div_value ~ group, data = T2_beta_df_chemical_combined)
m3<-aov(beta_div_value ~ group, data = T3_beta_df_chemical_combined)
summary(m3)
TukeyHSD(m3, "group")
boxplot(beta_div_value ~ group, data = T3_beta_df_chemical_combined)
beta_div__chemical_combined_T1_T2_T3 <- rbind(T1_beta_df_chemical_combined, T2_beta_df_chemical_combined, T3_beta_df_chemical_combined)
beta_anova(beta_div__chemical_combined_T1_T2_T3, Herbicide = "Chemical")
beta_anova(beta_div__chemical_combined_T1_T2_T3, Herbicide = "Hand")
beta_anova(beta_div__chemical_combined_T1_T2_T3, Herbicide = "Non-Treated")
treatment to control
plotDistances = function(p, m, s, d) {
# calc distances
wu = phyloseq::distance(p, m)
wu.m = melt(as.matrix(wu))
# remove self-comparisons
wu.m = wu.m %>%
filter(as.character(Var1) != as.character(Var2)) %>%
mutate_if(is.factor,as.character)
# get sample data (S4 error OK and expected)
sd = data.frame(sample_data(p)) %>%
select(s, d) %>%
mutate_if(is.factor,as.character)
sd$Herbicide <- factor(sd$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
# combined distances with sample data
colnames(sd) = c("Var1", "Type1")
wu.sd = left_join(wu.m, sd, by = "Var1")
colnames(sd) = c("Var2", "Type2")
wu.sd = left_join(wu.sd, sd, by = "Var2")
#remove this line to plot all comparisons.
wu.sd = wu.sd %>% filter(Type1 == "Hand" | Type1 == "Non-Treated")
# plot
ggplot(wu.sd, aes(x = Type2, y = value)) +
theme_bw() +
geom_point() +
geom_boxplot(aes(color = ifelse(Type1 == Type2, "red", "black"))) +
scale_color_identity() +
facet_wrap(~ Type1, scales = "free_x") +
theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
ggtitle(paste0("Distance Metric = ", m))
}
a<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T1"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
a <- a + ggtitle("Time 1 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T1_rare_allgroup_beta.pdf")
b<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T2"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
b <-b + ggtitle("Time 2 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T2_rare_allgroup_beta.pdf")
c<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T3"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
c<- c + ggtitle("Time 3 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T3_rare_allgroup_beta.pdf")
library(ggpubr)
ggarrange(a, b, c, ncol = 1)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_allgroup_beta.pdf", width = 7, height = 12)
Taxon abundance bar plot
#create super long color vector
col_vector <- c("#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941", "#006FA6", "#A30059",
"#FFDBE5", "#7A4900", "#0000A6", "#63FFAC", "#B79762", "#004D43", "#8FB0FF", "#997D87",
"#5A0007", "#809693", "#FEFFE6", "#1B4400", "#4FC601", "#3B5DFF", "#4A3B53", "#FF2F80",
"#61615A", "#BA0900", "#6B7900", "#00C2A0", "#FFAA92", "#FF90C9", "#B903AA", "#D16100",
"#DDEFFF", "#000035", "#7B4F4B", "#A1C299", "#300018", "#0AA6D8", "#013349", "#00846F",
"#372101", "#FFB500", "#C2FFED", "#A079BF", "#CC0744", "#C0B9B2", "#C2FF99", "#001E09",
"#00489C", "#6F0062", "#0CBD66", "#EEC3FF", "#456D75", "#B77B68", "#7A87A1", "#788D66",
"#885578", "#FAD09F", "#FF8A9A", "#D157A0", "#BEC459", "#456648", "#0086ED", "#886F4C",
"#34362D", "#B4A8BD", "#00A6AA", "#452C2C", "#636375", "#A3C8C9", "#FF913F", "#938A81",
"#575329", "#00FECF", "#B05B6F", "#8CD0FF", "#3B9700", "#04F757", "#C8A1A1", "#1E6E00",
"#7900D7", "#A77500", "#6367A9", "#A05837", "#6B002C", "#772600", "#D790FF", "#9B9700",
"#549E79", "#FFF69F", "#201625", "#72418F", "#BC23FF", "#99ADC0", "#3A2465", "#922329",
"#5B4534", "#FDE8DC", "#404E55", "#0089A3", "#CB7E98", "#A4E804", "#324E72", "#6A3A4C",
"#83AB58", "#001C1E", "#D1F7CE", "#004B28", "#C8D0F6", "#A3A489", "#806C66", "#222800",
"#BF5650", "#E83000", "#66796D", "#DA007C", "#FF1A59", "#8ADBB4", "#1E0200", "#5B4E51",
"#C895C5", "#320033", "#FF6832", "#66E1D3", "#CFCDAC", "#D0AC94", "#7ED379", "#012C58",
"#7A7BFF", "#D68E01", "#353339", "#78AFA1", "#FEB2C6", "#75797C", "#837393", "#943A4D",
"#B5F4FF", "#D2DCD5", "#9556BD", "#6A714A", "#001325", "#02525F", "#0AA3F7", "#E98176",
"#DBD5DD", "#5EBCD1", "#3D4F44", "#7E6405", "#02684E", "#962B75", "#8D8546", "#9695C5",
"#E773CE", "#D86A78", "#3E89BE", "#CA834E", "#518A87", "#5B113C", "#55813B", "#E704C4",
"#00005F", "#A97399", "#4B8160", "#59738A", "#FF5DA7", "#F7C9BF", "#643127", "#513A01",
"#6B94AA", "#51A058", "#A45B02", "#1D1702", "#E20027", "#E7AB63", "#4C6001", "#9C6966",
"#64547B", "#97979E", "#006A66", "#391406", "#F4D749", "#0045D2", "#006C31", "#DDB6D0",
"#7C6571", "#9FB2A4", "#00D891", "#15A08A", "#BC65E9", "#FFFFFE", "#C6DC99", "#203B3C",
"#671190", "#6B3A64", "#F5E1FF", "#FFA0F2", "#CCAA35", "#374527", "#8BB400", "#797868",
"#C6005A", "#3B000A", "#C86240", "#29607C", "#402334", "#7D5A44", "#CCB87C", "#B88183",
"#AA5199", "#B5D6C3", "#A38469", "#9F94F0", "#A74571", "#B894A6", "#71BB8C", "#00B433",
"#789EC9", "#6D80BA", "#953F00", "#5EFF03", "#E4FFFC", "#1BE177", "#BCB1E5", "#76912F",
"#003109", "#0060CD", "#D20096", "#895563", "#29201D", "#5B3213", "#A76F42", "#89412E",
"#1A3A2A", "#494B5A", "#A88C85", "#F4ABAA", "#A3F3AB", "#00C6C8", "#EA8B66", "#958A9F",
"#BDC9D2", "#9FA064", "#BE4700", "#658188", "#83A485", "#453C23", "#47675D", "#3A3F00",
"#061203", "#DFFB71", "#868E7E", "#98D058", "#6C8F7D", "#D7BFC2", "#3C3E6E", "#D83D66",
"#2F5D9B", "#6C5E46", "#D25B88", "#5B656C", "#00B57F", "#545C46", "#866097", "#365D25",
"#252F99", "#00CCFF", "#674E60", "#FC009C", "#92896B")
phylumGlommed <- tax_glom(ps_rare, "Phylum")
#t1
phylumGlommed_herb_t1 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T1"), group = "Herbicide")
phylumGlommed_herb_t1 <- transform_sample_counts(phylumGlommed_herb_t1, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t1)$Herbicide <- factor(sample_data(phylumGlommed_herb_t1)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t1, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_Taxon_barplot_t1.pdf")
#t2
phylumGlommed_herb_t2 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T2"), group = "Herbicide")
phylumGlommed_herb_t2 <- transform_sample_counts(phylumGlommed_herb_t2, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t2)$Herbicide <- factor(sample_data(phylumGlommed_herb_t2)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t2, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/ITS_Taxon_barplot_t2.pdf")
#t3
phylumGlommed_herb_t3 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T3"), group = "Herbicide")
phylumGlommed_herb_t3 <- transform_sample_counts(phylumGlommed_herb_t3, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t3)$Herbicide <- factor(sample_data(phylumGlommed_herb_t3)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t3, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/ITS_Taxon_barplot_t3.pdf")
Combined herbicide and time bar plot
sample_data(ps_rare_sub)$herb_time<-paste(sample_data(ps_rare_sub)$Herbicide, sample_data(ps_rare_sub)$Time, sep = "_")
ps_rare_for_barplot <- prune_taxa(taxa_sums(ps_rare_sub) > 50, ps_rare_sub)
plot_bar(ps_rare_for_barplot, x= "herb_time", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_BarPlot_Herbicide_Time.pdf", width = 20, height = 11)
Linear modeling of abundant taxa
mod_abund(test, IV = "Time")
Df Sum Sq Mean Sq F value Pr(>F)
matrix(data[, IV]) 2 30.58 15.292 3.725 0.0378 *
Residuals 26 106.73 4.105
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
[1] "Dicamba" "f.Coniochaetaceae"
Df Sum Sq Mean Sq F value Pr(>F)
matrix(data[, IV]) 2 864.4 432.2 7.708 0.00207 **
Residuals 29 1626.0 56.1
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
[1] "Handweeded" "f.Spizellomycetaceae"
Df Sum Sq Mean Sq F value Pr(>F)
matrix(data[, IV]) 2 805 402.7 3.351 0.0496 *
Residuals 28 3366 120.2
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
[1] "Non-Treated" "f.Spizellomycetaceae"
Df Sum Sq Mean Sq F value Pr(>F)
matrix(data[, IV]) 2 397.2 198.58 4.159 0.0255 *
Residuals 30 1432.4 47.75
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
[1] "Glyphosate" "f.Ceratobasidiaceae"
Df Sum Sq Mean Sq F value Pr(>F)
matrix(data[, IV]) 2 204.1 102.03 9.469 0.000649 ***
Residuals 30 323.3 10.78
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
[1] "Glyphosate" "f.Spizellomycetaceae"